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I .   INTRODUCTION 

A  program  is  underway  at  the  Naval  Postgraduate  School's 
(NPS)  Turbopropulsion  Laboratory  to  evaluate  the  wave  rotor 
concept.   The  wave  rotor,  operating  as  a  component  in  a  gas 
turbine  engine,  uses  unsteady  wave  propagation  in  tube-like 
passages  to  compress  incoming  air  before  it  goes  to  the  com- 
bustion chamber.   The  combustion  chamber  output  is  routed  back 
to  the  rotor  to  create  the  unsteady  waves .   Thus  the  rather 
simple  rotor,  with  "partial  admission"  inlet  and  outlet  ports, 
acts  both  as  a  compressor  and  as  a  turbine.   The  alternation 
of  hot  and  cool  gases  through  the  same  hardware  aids  cooling 
and  allows  higher  operating  (cycle)  temperatures  to  be  used. 
The  interaction  of  the  hot  and  cool  gases  within  the  passages 
of  the  wave  rotor  through  the  wave  patterns  that  are  created 
pose  a  difficult  problem.   The  work  underway  aims  to  develop 
and  validate  preliminary  design  and  performance  analysis  tools 

One  of  the  tools  needed  is  a  computer  program  which  can 
be  used  to  construct  wave  interactions  in  a  design  process 
and  then  accurately  predict  the  performance  of  the  design. 
Until  efficient  and  accurate  methods  of  solving  the  full 
Navier-Stokes  equations  for  unsteady  turbulent  flows  are 
developed,  the  design  program  will  be  based  initially  on  the 
solution  of  the  unsteady  Euler  equations  and  a  method  devised 
to  represent  losses.   Once  the  program  is  developed,  it  must 
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be  verified  against  experiment.   This  is  the  overall  goal  of 
the  present  program  in  which  a  wave  rotor  apparatus  has  been 
assembled  and  methods  of  measuring  the  unsteady  pressures  and 
temperatures  are  being  developed  concurrently  with  the  compu- 
tational effort. 

Three  different  approaches  to  the  solution  of  the  unsteady 
Euler  equations  were  examined  in  the  overall  program.   First, 
Eidelman  developed  a  tv/o-dimensional  code  based  on  the 
Godunov  method  of  solution  [Ref.  1]  and  applied  the  code 
to  examine  unsteady  wave  propagation  in  ducts  [Ref.  2]  and 
the  process  of  port  opening  to  wave  rotor  passages  [Ref.  3] .   A 
summary  of  Eidelman 's  work  is  given  in  Reference  4.   While 
the  method  is  conservative  and  does  not  require  the  introduc- 
tion of  artificial  viscosity,  the  extension  from  one  to  two 
dimensions  is  not  rigorous  when  shock  waves  are  present,  and 
computational  tim.es  with  the  present  code  are  quite  long. 
The  extension  to  include  viscous  effects  would  require  a  separate 
treatment  of  the  boundary  layer. 

Second,  Mathi.ir  developed  a  one-dimensional  code  based  on 
the  Random  Choice  Method  (RCM)  of  solution  [Ref.  5].   Somewhat 
similar  to  the  Godunov  method,  in  that  the  solution  is  based 
on  solving  the  Riemann  problem  within  each  grid  cell  at  each 
time  step,  the  RCM  approach  results  in  very  sharp  discontinui- 
ties which  can  be  tracked  easily.   This  is  particularly  useful 
in  constructing  wave  rotor  cycles,  in  which  the  position  of 
the  gas-gas  interface  is  equally  as  important  as  the  position 
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of  the  compression  and  expansion  waves.   The  code  is  there- 
fore valuable  in  the  preliminary  design  process,  to  examine 
suitable  port  arrangements,  and  the  gas  properties  at  the 
ports,  for  a  given  task.   Unfortunately,  an  extension  to  two- 
dimensional  and/or  viscous  flov;  can  not  be  made  rigorously. 

The  third  approach  was  followed  in  the  present  work.   The 
QAZID  method  for  compressible  inviscid  flow  computations 
developed  by  Verhoff  and  O'Neil   [Ref.  6]  was  implemented  to 
generate  a  one-dimensional  unsteady  Euler  code  with  the  goal 
of  evaluating  the  suitability  of  the  approach  for  wave  rotor 
applications.   Some  advantages  of  the  QAZID  method  were  recog- 
nized as  being  the  following: 

1.  The  method  is  based  on  the  use  of  characteristics. 
Such  methods  can  model  wave  propagation  accurately. 

2.  The  use  of  a  natural  streamline  coordinate  system 
eases  the  difficult  task  of  computing  with  two  and 
three  dimensional  grids. 

3.  The  equations  are  written  in  a  form  which  allows  a 
straightforward  extension  to  viscous  flows. 

4.  Codes  for  computing  internal,  steady,  two-dim.ensional 
flows  in  the  presence  of  shocks,  and  simple  internal 
viscous  flows,  have  been  generated  quickly  without 
significant  development  problems  [Ref.  6]. 

In  the  work  reported  herein,  a  one-dimensional  FORTRAN 
code  based  on  the  QAZID  method  was  developed  using  the  NFS 
IBM370-3033  computer  and  subsequently  exercised  on  the  shock 
tube  test  problem.   In  reporting  the  work,  first,  in  Section 
II  and  Appendix  A,  a  complete  account  is  given  of  the  deri- 
vation and  non-dimensionalization  of  the  governing  equations. 
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In  Section  III,  a  FORTRAN  code,  EULER-1,  is  described.   Its 
operation  on  the  NPS  computer,  and  the  listing  of  the  code, 
are  given  in  Appendix  B  and  Appendix  C,  respectively. 
Results  of  applying  the  code  to  the  shock  tube  test  problem 
are  given  in  Section  IV.   Difficulties  encountered  in  the 
implementation  of  the  method  and  additional  comments  are  given 
in  Section  V,  and  conclusions  are  given  in  Section  VI. 
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II .   THE  QAZID  METHOD 

A.  OVERVIEW 

The  QAZlD  method  uses  Riemann-like  variables  with  a  modi- 
fied entropy  term  in  their  definitions,,  to  express  the  Euler 
equations  in  a  natural  streamline  coordinate  system.   The 
resulting  partial  differential  equations  (PDE) ,  when  cast 
along  the  characteristic   trajectories  in  the  space-time 
domain,  reduce  to  a  system  of  ordinary  differential  equations 
(ODE)  which  may  be  solved  explicitly.   The  advantage  in  using 
the  modified  Riemann  variables  is  that  they  are  less  affected 
by  discontinuities  in  the  flow  than  are  the  standard  Riemann 
variables.   Since  the  equations  are  not  valid  across  discon- 
tinuities which  cause  irreversible  losses,-  such  discontinui- 
ties m.ust  be  located  and  treated  with  special  logic  in  the 
numerical  formulation. 

In  the  following  presentation,  the  reader's  familiarity 
with  the  method  of  characteristics  is  assumed. 

B.  DEVELOPMENT  OF  THE  EQUATIONS 

This  section  presents  the  governing  equations  and  outlines 
their  development.   A  rigorous  derivation  of  the  equations 
is  presented  in  Appendix  A= 

1.   The  Coordinate  System 

The  natural  streamline  coordinate  system  fs,n,m),  is 
shown  in  Fig.  A-1  relati"^^e  to  a  fixed  rectangular  cartesian 
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system  (x,y,z).   The  (s,n,m)  system  is  a  right-hand  orthogonal 
system  which  undergoes  curvilinear  translation  as  it  m.oves 
with  a  fluid  particle  along  a  streamline.   The  coordinate 
system  is  described  in  detail  in  Appendix  A. 
2 .   Variables 

The  modified  Riemann  variables,  or  "extended  Riemann 
variables"  [Ref.  6:p.  1],  are  defined  as 

Q   =   q  +  AS 

(1) 
R   =   q  -  AS 

where  q  is  the  velocity  magnitude,  A  is  the  speed  of  sound 
and  S  is  the  modified  entropy  defined  in  terms  of  pressure 
and  density  as 


s     =    7[7riyt2Y  -  £n(p/pM]  (2) 


The  modified  entropy  relation  is  the  result  of  defining  the 
entropy  change,-  dS ,  to  be  given  by 


1  '^Qr 

dS   =   -  -  — ^  (3 

Y   T 


where  dQ   is  the  heat  required  to  be  added  in  a  reversible 
process  between  the  same  end  states,  T  is  the  temperature  and 
Y  is  the  ratio  of  specific  heats. 


3 .  Conservation  of  Mass 

By  app.lying  the  continuity  equation  to  a  differential 
stream  tube  of  variable  cross  sectional  area,  in  natural 
coordinates,  and  ignoring  all  third  order  and  higher  deriva- 
tives, one  obtains 

|£  +  q  |£+  p  la  +  pq[||+  cos  e  |1]   =   0  (4) 

dt      dS      ds       dn  dm 

where  q  is  the  velocity,  p  is  the  density,  s  is  the  stream- 
wise  spatial  dimension,  and  9  and  cf)  are  the  flow  angles  as 
defined  in  Fig.  A-1  of  Appendix  A. 

4 .  Conservation  of  Momentum 

Applying  the  vector  form  of  the  momentum  conservation 
law  in  natural  coordinates,  the  equation  of  motion  for  inviscid 
flow  becomes 

:  r   9q  ^     9q  ^  9?. 
S^P  3t  ^  p^ai  ^  3¥^ 

n  ^  ^  9t    ^  ^  9s    9n 

+  i^[p  qcos  9  gf  +  p  q-cos  Q    J^   +    ^]       =       0 

5 .  The  Conservation  of  Energy 

In  the  absence  of  friction  and  heat  conduction,  and 
outside  of  discontinuities  (through  which  irreversible 
changes  consistent  with  mass,  momentum  and  energy  conservation 
are  permitted) ,  energy  conservation  is  equivalent  to  the 
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statement  that  the  entropy  of  a  fluid  particle  does  not  change, 
In  the  natural  coordinate  system,  therefore 


dt   =   8t  ^  ^  9¥   =   °  ^^^ 


Equation  (6)  will  not  be  valid  across  shock  waves  or  contact 
surfaces  between  gases  having  different  states. 
6 .   Transformation  to  a  Useful  Form 

Equations  (4)  and  (5)  are  first  expressed  in  terms 
of  the  primary  variables  q,  A,  S,  and  the  flow  angles.   Using 
the  definition  of  sound  speed  in  a  perfect  gas 


A^   =   yP/p  (7) 


Eq.  (4)  becomes 


at  "■  ^  ¥i  ^  ^Al¥  ^  ^/  2  fsH  ^  ^°^  ^  9lS^  =  0      ^'2) 

with  Eq.  (3)  ,  the  equation  of  state  for  a  perfect  gas,-  and 
the  first  law  of  thermodynamics,  Eq.  (5)  gives 


M+^iq  +  ^A_iA  +  A2^=o 

9t    ^  8s  ^  Y-1  8s    ^   3s      ^ 

le  +  ^  M   ^   _  A^  9  £n  P 

9t    ^  3s        Yq   8n  ^  ' 

11  +  a  ^   =   _     A^     9  g.n  P 
9t    ^  9s        Y  gcos  e   9m 
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Equations  (6),  (8)  and  (9)  constitute  the  system  of  P,D.E.'s 
which  describe  the  unsteady  isentropic  flow  of  a  perfect  gas 
in  natural  coordinates.   Since  third  order  terms  v^ere 
neglected  in  the  conservation  of  mass,  the  system  is  only 
accurate  to  second  order. 

The  system  of  P.D.E.'s  is  now  transformed  into  a  system 
of  OcD.E.'s  along  the  characteristic  trajectories.  In  general,- 
if  w  is  a  function  of  (s,t) ,  then 


dw   =   ^^  ds  +  ^^—  dt 
3s       dt 


or 


dw   _   ^fds_\  ,  ^ 
dt   ~   ^^dt^    Tt 


(10 


where  (ds/dt)  describes  the  "characteristic"  direction  in  the 
(s,t)  plane  along  which  Eq.  (10)  gives  the  rate  of  change  of 
the  parameter  w.   After  the  appropriate  algebra  and  the  intro- 
duction of  Eq .  (1),  the  final  system,  of  equations  becomes 


i  ^  <^-^'l§ 


9R 


dR 


9t    +     ^^-^)^ 


^  A  (s  -~)  [|T(q  --^  A  )] 

2  y-l        ds    ^       y-l 


.  izi.AsiM.cos  e||i 


=       +    ^A(S-^)[|^(q+^A)l 


+  IziqASlI^  *   COS    e   1^] 

Z  dn  dm. 


^ 


r 


(11 
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9S  ^    9S 
3t  ^  ^  9¥ 


96      99 
9t    ^  9s 


=   0 


A^  9  in   P 
yq   9n 


^  +   a  ^     =   - 


9  £n  P 


9t    ^  9s 


y   q  cos 9    9in 


The  characteristic  directions  in  the  space-time  plane  are 
clearly  q+A ,  q-A,  and  q. 

C.   SOLUTION  METHOD 

Equation  (11)  may  be  expressed  as 


f-  !^l  H  =  ^ 


(12) 


or  along  the  A  directions  as 


dw 
dt 


=   Z 


(13) 


where  the  vectors  are  defined  as 


^q' 


w   =    -s 


R 
S> 

9 


^ 


q+A 


q-A 


q 


Z  =  right  hand 
side  of 
Eq.  (11) 


22 


A  trajectory  in  the  space-time  domain  is  illustrated  in   Fig.  1, 


t+6t 


)■       6t 


[6w   =   v/g-Wg    ] 
i 


L Aw   =   w      -w^ ] 
s  •       A 


Figure  1.   Solution  Procedure  in  the  Space-Time  Domain 

which  shows  an  infinitesimal  interval  of  space  between  two 
"nodes"  of  the  computational  mesh. 

For  the  change  in  w  from  A  to  B  .along  the  trajectory 
with  slope  A, 


t+6t  _        _    _ 
/      Z  dt   =   ^B  "  "^A 


w_  -  w   J  +  [w   -  w, 
•  B    s  .        s  .    A" 

1  X 


=   6w  +  Aw 
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v/here  6w  denotes  the  change  due  to  time  at  a  fixed  location 
and  Aw  denotes  the  change  due  to  displacement  at  a  fixed 
time,  for  the  characteristic  trajectory.   The  essence  of  the 
solution  procedure  is  to  calculate  the  change  in  the  varia- 
bles at  each  spatial  node  during  the  differential  time 
interval,  so  that  the  step  to  the  next  time  level  can  be 
made.   In  other  words,  6w  must  be  calculated  at  each  node 
using 


_       _     t+6t 
6w   =   -Aw  +   /      Z  dt  (14 

t 


Since  all  the  information  at  time  level  t  is  known,  the 
position  of  A  can  be  determined  by  an  iterative  procedure 
based  on  X    and  Aw  can  be  calculated  by  interpolation.   The 
line  integral  can  be  evaluated  by  transforming  the  integral 
to  a  purely  spatial  integral  using  A  =  ds/dt.   Thus 


t+6t  B  -  ^i  - 

/      Z  dt   =    /   Y  ds   =       /    Y  ds 

t  A    '         s.+As 

1 


and  the  integral  can  be  evaluated  by  any  one  of  several 
numerical  methods.   Thus,  the  variables  at  each  node  can  be 
updated  in  time  using  Eq .  (14)  and  the  process  repeated. 
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D.   DISCONTINUITIES  IN  THE  FLOW 

There  are  several  types  of  discontinuities  that  must  be 
considered.   They  are  gradient  discontinuities,  contact 
discontinuities  and  shock  waves. 

Gradient  discontinuities  are  characteristic  of  the  head 
and  tail  of  rarefraction  waves,  the  collision  of  two  shocks  and 
the  interaction  of  a  shock  and  a  contact  discontinuity. 
Across  such  a  gradient,  the  derivatives  of  velocity,-  pressure 
and  sound  speed  are  discontinuous. 

A  contact  discontinuity  is  caused  by  the  interaction  of 
two  shocks  of  opposite  family  and  v/hen  a  shock  overtakes  a 
shock  of  the  same  family.   Entropy  and  sound  speed  are  discon- 
tinuous at  a  contact  discontinuity. 

Across  shock  waves,  velocity,  pressure,  density  and 
entropy  are  discontinuous  [Ref.  7]. 

Since  Eq.  (11)  is  not  valid  across  discontinuities,  addi- 
tional logic  is  required  in  the  numerical  procedure  to  model 
flows  which  contain  discontinuities.   The  method  presented  in 
Reference  6  for  making  a  correction  in  the  case  of  shock  waves 
will  give  good  accuracy  in  problems  where  the  solution  is 
converging  to  a  steady  state  condition  but  will  not  give 
accurate  results  during  the  transient  portion  of  such  problems 
or  for  problems  with  only  unsteady  solutions.   In  the  case  of 
applications  to  the  wave  rotor,  it  would  certainly  be  necessary 
to  know  the  locations  of  the  contact  surfaces  between  the 
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hot  and  cool  gases  as  well  as  the  locations  of  the  shock  and 
expansion  waves. 

The  methods  used  in  the  present  work  to  correct  for  dis- 
continuities are  those  of  Moretti  [Ref.  7].   The  present 
discussion  will  be  limited  to  the  treatm.ent  of  shock  waves. 
The  method  makes  use  of  the  analytical  relationship  between 
the  change  in  the  Riemann  variables  across  a  shock  and  the 
incoming  Mach  number  relative  to  the  shock  wave  (W)  which  is 
illustrated  in  Figs.  2  and  3.   The  relation  is  used  to 
determine  the  shock  speed  and  to  transform  the  problem  to 
a  steady  case  which  can  be  handled  using  normal  shock  rela- 
tions.  For  the  situation  depicted  in  Fig.  2  of  a  shock 
propagating  to  the  right  with  velocity  V   into  air  with 
velocity  q„,  and  with  the  high  pressure  side  to  the  left, 
where  A  denotes  the  left  side  of  the  shock  and  B  the  right, 


If  the  pressure  and  density  are  non-dimensionalized  by  the 
values  on  the  low  pressure  side  of  the  shock,  the  change  in 
the  extended  Riem.ann  variable  Q  is  given  by 


^A  '^B   ^   2(W^-1)     2   \   ,. 
Ag       (Y  +  1)W   ^  Y-l'^  -^ 


A  r   1        r  r  2y  ,,2   v-l.  ,  .  (y-1)W^  +  2/. 
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Figure    2.       Shock  Wave  with   High   Pressure    to    the    Left 
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Figure    3.       Shock   Wave   with   High   Pressure    to    the    Right 
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where 


^  =  TY^f2'^-^"i^^"'H^w2-i))''' 


This  exact  relationship  can  be  approximated  by  the  polynomial 
V^B       ,  2 


Ab 


=   -2.7574  +  3.1573W  -  0.2863W  (17) 


as  illustrated  in  Fig.  4  over  a  shock  strength  range  of  1.0 
to  4 . 0 .   It  should  be  noted  that  if  the  high  pressure  side 
were  to  the  right,  as  in  Fig.  3,  Eq.  (15)  would  become 


A      ^A    s 
Vv   = 


A         A 


and  the  left  side  of  Eq.  (16)  would  become 


The  procedure  then,  is  to  measure  the  change  in  Q  across  an 
interval  where  a  shock  is  known  to  exist.   This  value,  call 
it  AQ  ,  is  used  in  Eq.  (17)  to  get  an  approximation  for  the 
corresponding  value  of  W  and  then  the  exact  value  of  AQ,  call 
it  AQ  ,  is  calculated  using  Eq.  (16).   If  the  exact  value  of 
AQ  is  not  equal  to  the  value  measured  across  the  shock,  a 
new  value,  AQ,  is  calculated  according  to 


AQ.  ^   =   AQ.  +  (AQ   -  AQ^) 
^1+1       ^1      ^m     E 


and  entered  into  Eq.  (17)  to  obtain  another  value  of  W.   When 
the  exact  value  of  AQ  calculated  from  Eq.  (16)  equals  the 
measured  value  of  AQ,  W  is  known  to  be  correct  and  the  normal 
shock  relations  can  be  used  to  calculate  the  values  at  node  A. 
Also,  V   can  be  calculated  from  Eq .  (15)  and  used  to  track 
the  shock  during  the  next  time  interval. 
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Ill,   DESCRIPTION  OF  FORTRAN  PROGRAM  "EULER-1" 

A.  MACHINE  AND  LANGUAGE 

The  EULER-1  program  is  written  in  VS  FORTRAN  and  runs  on 
an  IBM  30  33,  System  37  0  computer,  however  the  code  is  simple 
and  small  enough  to  enter  and  run  on  a  mini  or  micro  computer 
in  a  Basic  language  if  one  is  willing  to  accept  significantly 
longer  run  times.   Table  1  at  the  end  of  this  section  is  a 
summary  of  the  editable  parameters  and  their  effect  on  the 
program. 

B.  CONVENTIONS  AND  BASIC  STRUCTURE 

EULER-1  is  a  first-order  one-dimensional  code  using  an 
evenly  spaced  numerical  grid.   All  values  are  double  precision 
except  those  used  in  the  graphics  routines  which  are  rounded 
to  single  precision. 

Each  subroutine  has  its  own  variables  space.   In  other 
words,  variables  are  not  shared  in  common  throughout  the 
program  but  must  be  passed  to  the  subroutine  being  called 
by  the  calling  routine.   In  all  cases,  however,  the  variables 
have  the  same  name  in  both  the  called  and  the  calling  routine 
so  there  should  be  no  confusion.   This  was  done  so  that 
arrays  could  be  dimensioned  at  execution  tim.e. 

The  program,  depicted  in  Fig.  5,  is  structured  aroi.md  a 
main  routine  which  serves  as  a  user  input  area,  sets  up  the 
problem  and  calls  five  subroutines  to  solve  the  problem  and 
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Figure    5,       EULER-1    Flowchart 
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Figure    5.        (Continued 
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output  the  solution.   The  five  subroutines  are  "TIME,"  "TRAK," 
"SWEEP,"  "JUMP"  and  then  one  of  several  output  routines  depend- 
ing on  user  desires ,   Output  options  include  two  types  of 
graphical  displays,  "PLOT"  and  "EXACT,"  and  one  tabular 
listing  routine,  "LIST."   When  "PLOT"  is  selected,  "BORDER" 
is  automatically  called  to  set  up  the  plotting  area.   There  are 
two  other  subroutines,  "RSHOCK"  and  "LSHOCK,"  which  are  called 
by  "SV7EEP"  as  needed. 

C.   SUBROUTINE  DESCRIPTIONS 

In  general,  each  subroutine  begins  with  a  heading  followed 
by  variable  definitions  v/here  appropriate,  variable  declara- 
tions and  array  dimensioning.   Most  variables  are  defined  in 
the  "MAIN"  routine  and  only  those  variables  which  were  not 
are  defined  in  the  subroutines  where  they  are  used. 
1.   The  "MAIN"  Routine 

This  routine  forms  the  main  structure  of  the  program 
and  includes  a  heading,  an  extensive  list  of  variable  defini- 
tions, a  user  input  area,  the  necessary  statements  to  make 
the  initial  value  assignments  to  variables,  and  the  call 
statements  for  the  ^^arious  subroutines. 

The  input  area  is  those  lines  of  the  program  (145-175) 
where  the  user  edits  the  program  to  establish  the  initial 
conditions,  mesh  size,  term.ination  criteria  and  output  options. 

The  initial  conditions  which  m.ay  be  modified  are  the 
pressure,  temperature  and  density  ratios  across  the  diaphragm, 
the  initial  velocity  of  the  fluid  on  each  side  of  the 
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diaphragm,  and  the  value  of  gaipina  (which  must  be  the  same  for 
both  sides) . 

All  velocities  are  non-dimensionalized  by  the  initial 
sound  speed  on  the  low  pressure  side  of  the  diaphragm,  and 
pressures  and  densities  are  non-dimensionalized  by  their 
initial  values  on  the  low  pressure  side  of  the  diaphragm. 

The  mesh  size  may  be  set  at  any  odd  number  and  the 
arrays  in  the  user  input  area  must  be  dimensioned  as  such. 

The  criteria  for  program  termination  is  the  number  of 
time  steps  com.puted,-  which  the  user  selects. 

The  output  options  are  determ.ined  by  the  value  of  the 
variable  GRAPHS.   A  value  of  zero  results  in  a  call  to  "LIST" 
which  writes  to  file  9  on  the' -user ' s  permanent  disk,  a  tabular 
listing  of  the  variable  arrays  and  discontinuity  locations. 
A  value  of  1  results  in  a  call  to  "BORDER"  and  "PLOT"  which 
create  a  plot  of  the  pressure,  density,  velocity  and  entropy 
distributions  as  in  Fig.  7.   A  ^^-alue  of  2  results  in  a  call 
to  "EXACT"  which  creates  a  plot  of  the  density  distribution 
compared  to  the  exact  solution  as  in  Fig.  8.   The  exact  values 
to  be  plotted  must  be  entered  in  the  user  input  area,  they  are 
not  calculated  by  the  program.   A  more  detailed  description 
of  the  various  outputs  is  found  under  the  appropriate  subrou- 
tine description.   The  frequency  with  which  output  is  created 
is  controlled  by  the  variable  SKIP  which  is  the  number  of  tim.e 
steps  between  calls  to  output  routines. 
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2.  The  "TIME"  Routine 

The  maximum  allowable  time  step  is  governed  by  the 
CFL  condition.   Simply  stated,  this  means  that  the  time  step 
must  be  small  enough  so  that  the  characteristic  trajectories 
rem.ain  within  one  spatial  interval  during  the  time  interval. 
The  minimum  time  step  is  calculated  by  computing 

Delt   =   H/ABS [Q+A] 

at  every  node,  where  H  is  the  spatial  interval,  and  selecting 
the  minimum  value  of  Delt. 

3.  The  "TRAK"  Routine 

Shock  locations  at  the  unknown  time  level  are  deter- 
m.ined  by  computing  the  shock  speed  at  the  known  time  level, 
as  outlined  in  Section  2,  and  multiplying  by  the  time  inter- 
val computed  by  "TIME."   The  time  step  is  reduced,  if  neces- 
sary, to  limit  shock  travel  to  one  spatial  interval.   The  node 
imm.ediately  to  the  right  of  the  shock  (upstream)  is  flagged 
for  later  use  by  "SWEEP"  and  "JUMP." 

In  calculating  the  shock  speed,  the  assumption  is  made 
that  the  conditions  immediately  adjacent  to  the  shock  are  the 
sam.e  as  the  conditions  at  the  nodes  to  the  left  and  right  of 
the  shock. 

4.  The  "SWEEP"  Routine 

The  "SWEEP"  routine  makes  the  necessary  calculations 
to  solve  Eq.  (14) .   This  involves  interpolation  at  the  known 
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time  level  for  the  values  of  Aw,  calculation  of  z,  and 
integration  of  z  for  each  node. 

To  calculate  Aw,  the  assumption  is  made  that  the 
characteristic  trajectories  are  straight  lines.   An  initial 
guess  of  the  characteristic  slope  (A)  is  made  and  by  forcing 
the  characteristic  to  pass  through  the  node  at  point  B  in 
Figc  1,  As  =  A  xdelt.   The  values  of  q  and  A  at  point  A  are 
found  by  interpolation  and  the  slope  of  the  characteristic 
is  then  calculated  using  the  values  of  q  and  A  at  point  A, 
This  slope  is  compared  with  the  initial  guess  and  if  the  tv;o 
are  not  in  agreement,  the  procedure  is  iterated  until  they 
converge.   Once  As  is  accurately  determined,  the  values  of  Aw 
can  be  calculated  by  interpolation ,   Linear  interpolation  is 
used  in  the  present  version  of  EULER-1.   This  procedure  is 
carried  out  for  each  characteristic  at  each  node. 


t+6t 


t  >f 
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Figure  6.   "RSHOCK"  Case 

A  deviation  from  the  above  procedure  is  necessary  in 
the  case  where  a  shock  exists  in  the  spatial  interval  to  the 
left  or  right  of  the  current  node.   When  a  shock  exists  to 
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the  rights  and  the  flow  is  subsonic,  as  in  Fig.  6,  care  must 
be  taken  not  to  interpolate  for  point  3  based  on  the  slope 
between  points  T+1  and  I,  which  would  not  be  accurate  due 
to  the  discontinuity.   In  such  a  case  "RSHOCK"  is  called  and 
the  interpolation  is  based  on  the  slope  between  point  T  and 
I-l.   Similarly,  when  the  shock  is  to  the  left  of  the  current 
node,  "LSHOCK"  is  called  and  the  interpolation  is  based  on 
the  slope  between  I+l  and  I. 

The  calculation  of  z  includes  the  calculation  of 
spatial  derivatives  of  q  and  A  for  characteristics  1  and  3 
of  Fig.  6.   In  general,  the  derivatives  associated  with 
characteristic  3  are  forward  differenced  and  those  associated 
with  characteristic  1  are  backward  differenced  in  keeping 
with  the  principle  of  domain  of  dependence.   If  the  flow  is 
supersonic  or  if  a  shock  exists  to  the  right  of  the  current 
node,  all  derivatives  are  backward  differenced.   If  a  shock 
exists  to  the  left  of  the  current  node,  all  derivatives  are 
forward  differenced.   Once  the  derivatives  are  known,  z  is 
calculated  from  Eq .  (11)  using  the  average  value  of  A  between 
points  1  and  I  or  3  and  I  as  appropriate.   Since  the  deriva- 
tives are  linear,  this  results  in  an  average  value  of  z 
over  the  same  interval. 

The  integration  of  z  is  transform.ed  from  a  time  to  a 
spatial  integration  as  described  in  Section  II  and  the 
trapezoidal  rule  is  used  to  carry  out  the  integration.   In 
EULER-1,  this  has  been  done  in  one  step  using  the  average 
value  of  z  described  above. 
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Equation  (14)  is  then  solved  by  addition  of  Aw  and 
the  results  of  the  integration,  and  the  resulting  6w  is 
stored  until  all  nodes  are  calculated.   After  all  nodes  have 
been  calculated,  the  variable  arrays  (w)  are  updated  for  the 
next  time  interval. 

5.  The  "LSHOCK"  Routine 

As  discussed  above,  the  "LSHOCK"  routine  modifies 
the  basic  EULER-1  procedure  at  an  interior  node  when  a  shock 
exists  to  the  left  of  the  node.   Interpolation  of  quantities 
to  the  left  of  the  node  are  computed  based  on  the  derivatives 
of  the  quantities  to  the  right  of  the  node.   Although  this 
assumes  that  the  derivatives  do  not  change  between  adjacent 
spatial  intervals,  it  is  necessary  to  avoid  taking  derivatives 
across  discontinuities  in  the  flow. 

6.  The  "RSHOCK"  Routine 

Similar  to  "LSHOCK,"  "RSHOCK"  bases  interpolations 
of  quantities  to  the  right  of  the  node  on  the  derivatives  of 
the  quantities  +-o  the  left  of  the  node  when  a  shock  exists 
to  the  right  of  the  node. 

7.  The  "JUMP"  Routine 

The  "JUMP"  routine  is  used  to  calculate  the  conditions 
downstream  of  a  stationary  shock  as  described  in  Section  II. 
If  a  shock  is  known  to  have  crossed  a  spatial  node  during  a 
time  interval,  which  is  known  once  "TRAK"  has  been  called  for 
that  time  interval,  the  entire  "SWEEP"  sequence  is  skipped  for 
the  node  which  was  crossed  by  the  shock  and  the  conditions  at 
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that  node  are  determined  using  the  normal  shock  relations.   Note 
that  the  node  in  question  is  the  node  downstream  of  the  sta- 
tionary shock.   As  in  "TRACK,"  it  is  assumed  that  the  condi- 
tions at  the  nodes  upstream  and  downstream  of  the  shock  are 
the  same  as  the  conditions  immediately  adjacent  to  the  shock. 
8 .   The  Output  Routines 

There  are  four  subroutines  which  produce  various  types 
of  output  as  the  user  desires.   "BORDER"  and  "PLOT"  produce 
a  graphical  presentation  of  the  pressure,  density,  velocity 
and  entropy  distributions  at  selected  time  levels  as  seen 
for  example,  in  Fig.  7.   "EXACT"  produces  a  plot  of  the  den- 
sity distribution  at  one  selected  time  interval  and  compares 
it  with  the  exact  solution  as  shown  in  Fig.  8.   At  selected 
time  levels,  "LIST"  produces  a  tabular  listing  of  the  Riemiann 
variables,  modified  entropy,  pressure,  density  and  ^^^elocity 
distributions,  elapsed  time,  shock  speed  and  discontinuity 
locations.   The  listing  is  written  to  the  user's  permanent 
storage  disk. 

When  the  value  of  GRAPHS  is  set  equal  to  one  in  the 
"MAIN"  routine,  "BORDER"  is  called  once  to  set  up  the  plot 
axis,  labels,  and  headings.   "PLOT"  is  called  every  SKIP 
time  steps  to  draw  the  four  distrd.bution  curves. 

When  the  value  of  GRAPHS  is  set  equal  to  two,  "EXACT" 
is  called  every  SKIP  tim.e  steps  to  plot  the  density  distri- 
bution compared  with  the  exact  solution  computed  at  six 
points  of  interest.   The  points  are  the  two  end  points,  which 
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are  simply  the  initial  conditions,  the  head  and  tai].  of  the 
rarefraction  wave,  the  point  just  left  of  the  contact  surface 
and  the  point  just  left  of  the  shock.   The  spatial  locations 
of  these  points  are  computed  by  "EXACT"  based  on  the  elapsed 
tim.e  and  the  known  values  of  the  wave  velocities  entered  in 
the  "MAIN"  routine.   The  exact  values  of  the  densities  at  these 
points  must  also  be  entered  in  the  "MAIN"  routine. 

When  the  value  of  GRAPHS  is  set  equal  to  zero,  the 
tabular  listing  as  described  above  is  sent  to  the  user's 
disk.   No  graphical  output  is  created. 

D.   CAPABILITIES  AND  LIMITATIONS 

The  present  version  of  EULER-1  is  set  up  to  solve  a  shock 
tube  problem  with  a  single  centered  diaphragm.   Boundary 
conditions  for  the  ends  of  the  tube  have  not  been  incorporated 
so  the  problem  must  be  stopped  before  the  waves  reach  the  end 
of  the  tube. 

The  left  side  of  the  diaphragm  is  the  high  pressure  side. 
The  program  wil.l  not  run  with  the  right  side  as  the  high 
pressure  side  without  changes  to  some  of  the  shock  correction 
logic. 

The  user  may  select  any  odd  number  of  grid  points  limited 
only  by  the  am.ount  of  memory  available. 

The  program  tracks  shock  waves  and  makes  shock  jump  calcu- 
lations at  the  appropriate  locations  but  the  program  can  also 
run  without  the  shock  tracking  feature  and  jump  calculations 
if  so  desired  with  some  smearing  of  the  shock  discontinuity 
and  complete  loss  of  entropy  change  across  the  shock. 
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TABLE  1 
LIST  OF  EDITABLE  PARAMETERS 

Parameter  Function 

N  Number  of  grid  points 

M  Controls  which  discontinuities  are 

tracked 

1  =  Shocks  only 

2  =  Shocks  and  contact  surfaces 

GRAPHS  Controls  the  form  of  the  output 

0  =  Tabular  listing 

1  =  Pressure,  density,  velocity, 

entropy  plot 

2  =  Exact  solution  comparison  for 

density 

SKIP  Number  of  time  steps  betv/een  output 

calls 

JSTOP  Number  of  time  steps  calculated 

TRI  Initial  temperature  ratio 

PRI  Initial  pressure  ratio 

DRI  Initial  density  ratio 

QLI  Initial  velocity  left  of  the  diaphragm 

QRI  Initial  velocity  right  of  the  diaphragm 

G  Ratio  of  specific  heats 
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IV.   TEST  RESULTS 

EULER-1  was  tested  on  the  Riemann  shock  tube  problem. 
The  initial  conditions  and  the  courseness  of  the  computa- 
tional mesh  were  varied.   Additionally,  one  run  was  made 
without  the  shock  tracking  and  correction  features. 

Figures  7  and  8  show  the  results  for  initial  pressure 
and  density  ratios  of  5,  uniform  temperature  and  with  the  air 
initially  at  rest,  using  a  grid  of  101  points.   The  exact, 
analytically  predicted  conditions  are  also  shown  in  Fig.  8. 
It  can  be  seen  that  all  wave  velocities  are  correctly  computed 
and  the  shock  is  defined  within  one  interval .   The  rarefrac- 
tion  waves  are  slightly  smeared  and  the  contact  surface  is 
greatly  smeared.   A  slight  transient  instability  to  the  right 
of  the  diaphragm  is  noticeable,  particularly  in  the  plots 
of  pressure  and  velocity. 

Figures  9  and  10  show  results  for  initial  pressure  and 
density  ratios  of  1.3.   The  shock  is  still  well  defined  in 
the  correct  interval  and  the  contact  surface  is  not  as  badly 
smeared  as  in  Fig.  8.   There  is  a  loss  of  accuracy  however, 
with  respect  to  the  tail  of  the  rarefraction  wave. 

Figure  11  shows  the  results  of  the  original  problem  with 
a  grid  of  51  points.   As  might  be  expected  with  a  coarser 
mesh,  the  discontinuities  are  less  sharply  defined  although 
it  is  clear  that  the  wave  ^^elociMes  are  accurately  computed. 


42 


Figure  12  shows  the  results  for  the  original  problem  with 
a  grid  of  901  points.   It  can  be  seen  that  the  EULER-1 
solution  closely  approximates  the  exact  solution. 

Figures  13  and  14  show  the  result  of  running  the  original 
problem  without  the  shock  tracking  and  correcting  features. 
Note  that  the  shock  position  is  incorrectly  computed  and  it 
is  smeared  slightly.   Also  note  in  Fig.  13,  the  lack  of  any 
change  in  modified  entropy  across  the  shock. 

Run  time  for  the  101  point  mesh  is  0.00049  seconds  per 
node  per  time  step.   The  results  for  Fig.  7  took  50  time 
steps  for  a  total  time  of  2.46  seconds.   For  the  901  point 
mesh,  the  run  time  decreased  to  0.00023  seconds  per  node  per 
time  step  due  to  the  less  frequent  use  of  the  shock  tracking 
and  correction  routines  per  node  calculated.   The  results  for 
Fig.  12  took  470  time  steps  for  a  total  time  of  71.85  seconds 
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SHOCK   TUBE   RESULTS 
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Figure  7.   EULER-1  Results  for  Pressure  Ratio  5  with 
Shock  Tracking 
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SHOCK   TUBE   RESULTS 
FIRST   ORDER  N   =    101 

DENSITY   RATIO   =    1.3        TEMP  RRTIO 
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Figure  9.   EULEP-1  Results  for  Pressure  Ratio  1,3  with 
Shock  Trackinq 
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SHOCK   TUBE   RESULTS 
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Figure  13.   EULER-1  Results  for  Pressure  Ratio  5 
without  Shock  Tracking 
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V.   DISCUSSION 

A.  RESULTS 

The  transient  instability  noted  in  the  results  was  not 
investigated  since  the  amplitude  was  small  and  the  effect 
would  be  of  no  consequence  in  most  applications. 

At  low  pressure  ratios,  the  velocity  of  propagation  of 
the  tail  of  the  rarefraction  wave  approached  sound  speed 
and,  in  these  cases,  the  EULER-1  code  did  not  accurately 
compute  its  location.   The  reason  for  this  was  not  investi- 
gated here;  however,  due  to  the  low  pressure  ratios  at  which 
wave  rotors  can  operate  i   this  is  of  interest  and  should  be 
investigated  in  future  work. 

B.  SPECIAL  CONSIDERATIONS 
1 ,   Modified  Entropy 

In  the  QAZID  method,  and  subsequently,  in  the  EULER-l 
code,  the  variable  S,  which  here  has  been  referred  to  as  the 
"modified  entropy"  and  which  is  referred  to  in  Reference  6  as 
simply  the  "entropy,"  does  not  behave  thermodynamically  as 
one  would  expect  entropy  to  behave.   As  a  fluid  crosses  a 
shock  wave,  the  modified  entropy  of  the  fluid  decreases. 
This  is  the  result  of  the  definition  of  Eq .  (3) ,  repeated  here, 
that 

1  ^Qr 

dS   =   -  -  — 
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The  use  of  the  word  "entropy"  and  symbol  "S"  for  this  varia- 
ble is  confusing  since  the  sign  of  the  change  in  "S"  is  in 
direct  conflict  with  well-established  conventions  in  thermo- 
dynamics.  For  example,  Corollary  6  to  the  2    Law  of 
Thermodynamics  [Ref.  8:p,  88]  would  have  to  be  changed  to 
read,  "The  'modified  entropy'  of  an  isolated  system  decreases 
or  in  the  limit  remains  constant."   At  a  minimum,  perhaps,  the 
symbol  §  to  denote  the  "modified  (negative  and  scaled)  entropy" 
would  be  helpful.   In  practice  of  course,  the  change  in  the 
conventionally  defined  entropy  can  be  recovered  from  the  change 
in  the  modified  entropy,  simply  by  multiplying  by  (-7). 
2 .   Moretti's  Methods 

Incorporating  Moretti's  methods  of  handling  discon- 
tinuities [Ref.  7]  into  the  EULER-1  code  required  special 
care. 

There  is  a  fundamental  difference  in  procedure  in 
that  the  EULER-1  code  is  based  on  the  high  pressure  side 
being  on  the  left  whereas  Moretti's  formulations  are  all 
based  on  the  high  pressure  side  being  on  the  right.   Although 
essentially  a  matter  of  bookkeeping,  it  is  an  area  of  poten- 
tial confusion. 

Another  difference  is  that  the  Riemann  variables  used 
by  Morreti  are  not  the  Riemann  variables  used  in  the  QAZID 
method.   This  can  lead  to  difficulty.   In  fact,  in  the 
application  of  Moretti's  shock  tracking  scheme  to  the  EULER-1 
code,  the  pressures  and  densities  have  to  be  non-dimensional ized 
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by  their  values  on  the  low  pressure  side  of  the  shock  in  order 
to  obtain  Eq.  (14),  and  for  the  procedure  to  work. 

C.   SUITABILITY  FOR  WAVE  ROTOR  APPLICATIONS 

Further  work  is  necessary  before  the  suitability  of  the 
QAZID  method  for  wave  rotor  applications  can  be  fully  evalu- 
ated.  The  following  areas  need  to  be  addressed. 

1 .  Boundary  Conditions 

In  the  EULER-1  code,  boundary  conditions  at  the  ends 
of  the  passage  have  not  been  incorporated.   The  code  calculates 
the  changes  of  the  interior  nodes  only,  skipping  the  two  end 
nodes.   In  particular,  solid  wall  boundary  conditions  and 
open-end  boundary  conditions  must  be  incorporated.   No  particu- 
lar problems  are  expected  in  the  boundary  conditions  themselves 
However,  additional  logic  may  be  necessary  in  the  handling  of 
discontinuities,  such  as  the  reversal  of  conditions  to  the 
left  and  right  of  a  discontinuity  after  reflection  from  a 
boundary.   The  additional  logic,  which  may  be  considerable, 
is  warranted  by  the  simplicity  of  the  QAZID  method. 

2.  Contact  Discontinuity 

No  special  attention  is  given  to  the  contact  discon- 
tinuity in  the  EULER-1  code  and  hence  the  discontinuity  is 
smeared.   In  a  wave  rotor  application,  this  discontinuity 
will  have  to  be  sharp.   Moretti's  methods  again  appear  to  he 
applicable  here  and  the  additional  logic  needs  only  to  be 
formulated  and  incorporated  into  the  EULER-1  code.. 
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3 .  Quasi-One-Dimensional  Modeling 

EULER-1  is  a  one-dimensional  code.   The  solution  of 
flows  in  passages  of  varying  area  may  be  desirable  for  some 
wave  rotor  applications.   Such  problems  can  be  solved  in  a 
quasi-one-dimensional  manner  by  the  addition  of  the  appro- 
priate area  variation  term  to  the  equations,  and  the  need  to 
solve  the  fully  two-dimensional  equations  is  avoided.   The 
area  variation  must  be  incorporated  into  the  EULER-1  code  and 
the  code  tested  on  quasi-one-dimensional  problems. 

4 .  Other  Potential  Extensions 

Another  capability  which  can  be  incorporated  is  the 
ability  to  handle  two  gases  with  different  specific  heat 
ratios.   Eventually,  the  code  should  be  extended  also  to  a 
second  order,  one-dimensional  version  and  then  to  a  first 
order,  two-dimensional  version. 
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VI.   CONCLUSIONS 

The  development  of  the  EULER  equations  in  the  natural 
streamline  coordinates  using  extended  Riemann  variables  was 
reviewed  in  detail.   The  QAZID  solution  method,  with  the 
addition  of  shock  tracking  features,  was  implemented  in  a 
first  order,  one-dimensional  FORTRAN  code  with  the  intent  of 
evaluating  the  method's  suitability  for  wave  rotor  applica- 
tions.  The  code  (EULER-1)  was  tested  on  the  shock  tube 
problem,  with  good  results.   The  incorporation  of  boundary 
conditions,  an  improvement  in  the  contact  discontinuity 
definition  and  the  addition  of  an  area  variation  term  for 
quasi-one-dimensional  modeling  are  considered  to  be  neces- 
sary before  the  QAZID  method  can  be  accepted  as  being  suitable 
for  wave  rotor  applications.   The  additional  logic  required 
for  these  extensions  may  be  considerable  but  the  development 
is  warranted  by  the  overall  simplicity  of  the  QAZID  method, 
and  its  straightforward  extension  to  viscous,  multi-dimensional 
flow  modeling. 
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APPENDIX  A 

DERIVATION  OF  THE  GOVERNING  EQUATIONS  FOR 
THREE-DIMENSIONAL  INVISCID  FLOW  (WITH  AREA  CHANGE) 

The  governing  equations  of  a  3-D,  inviscid,  compressible, 
unsteady  flow  are  derived  here  in  a  natural  streamline 
coordinate  system.   The  equations  are  then  recast  along 
characteristic  trajectories  in  the  (s,t)  plane  and  expressed 
in  extended  Riemann  variables,  reducing  the  system  of  par- 
tial differential  equations  to  a  system  of  ordinary  differ- 
ential equations. 

A.   DEFINITION  OF  VARIABLES 

A        Cross-sectional  area  of  a  differential  stream 
tube 

A       Speed  of  sound 

C        Specific  heat  at  constant  pressure 

C        Specific  heat  at  constant  volume 

V 

da       Normal  vector  for  differential  area  of  control 
surface 

e        Specific  internal  energy 

i        Unit  vectors  in  the  s ,-  n,  and  m  directions 
s ,  n,m 

P  Static  pressure 

Q  Modified  Riem.ann  variable 

Q  Re^^'ersibl.e  heat  transferred 

q  Velocity  magnitude 

r  Position  vector 
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R  Modified  Riemann  variable 

R  Gas  constant 

S  Modified  entropy 

T  Static  temperature 

t  Time 

u  Velocity  relative  to  a  standing  shock  wave 

V  Specific  volume 

vol  Differential  element  of  volume 

V  Velocity  vector 

W        Incom^ing  Mach  number  relative  to  a  standing 
shock  wave 

p       Density 

V  Ratio  of  specific  heats 

Q  ,  (p  Flow  angles  with  respect  to  reference  coordinate 

planes 

V  Vector  operator 

B,   THE  COORDINATE  SYSTEM 

The  natural  streamline  coordinate  system  (s,n,m)  is 
shown  with  respect  to  a  fixed  rectangular  cartesian  system 
(x,y,z)  in  Fig.  Al .   The  system  is  a  curvilinearly  translating, 
right  hand  orthogonal  system  which  translates  with  a  fluid 
particle  along  a  stream.line  such  that  the  's'  coordinate  is 
measured  in  the  direction  of  the  flow.   The  'n'  coordinate 
direction  always  lies  in  the  plane  defined  by  the  's'  coor- 
dinate direction  and  the  fixed  'y'  axis.   The  'm'  direction  is 
norma]  to  the  (s,n)  coordinate  surface.   The  flow  angles, - 
6  and  (|) ,  are  defined  as  shown  in  Fig.  Al . 


STREAMLINE 


>  X 


Figure  Al .   The  Natural  Coordinate  System 

C.   VECTOR  OPERATORS  IN  THE  NATURAL  COORDINATE  SYSTEM 
1.   ^(  ) 

If  dr  is  the  position  vector  of  a  fluid  particle 


dr  •  V(  )   =   d(  ) 


:ai 


In   natural   coordinates,  Eq,  (Al)  becomes 


ri   ds    +    i    dn    +    i    dm]-[(    )i      +    (    )i      +    (    )i    ] 
■s  n  m'-s  n  m 


'    >^s   +  ILlan   .  il^m 


ds  9n 


9m 


By    inspection. 
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9s    s    9n    n    3m    m 


2.   V-V(  ) 


Since 


V   =   q  i3 


using  Eq.  (A2) 


V-V(  )   =   q  1^  (A3 


3.   V-v 


From  Eq.  (A2)  with  V  =  qi 


V-V  =   If  (A4) 


4.    (V-V)V 


From  Eq.  (A3)  with  V  =  qi, 


(V.V)V   =   q  ^   =   q  -^^   =   q(Ji^  +  q  ^1  (AS: 


The  change  in  the  unit  vector  i   v/ith  respect  to  s  is 
derived  below  and  is  illustrated  in  Figs.  A2  through  A6 .   (A 
three-dimensional  model  is  very  helpfu],  in  visualizing  the 
vector  geometry.) 

Figure  A2  illustrates  the  natural  coordinate  system 
translating  along  a  streamline  from  point  A  to  point  B  in 
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m   rq 


■^  X 


< 


X 


Figure  A2 .   Unit  Vectors  Moving  Along  A  Streamline 


Y 


Figure  A3.   Superposition  of  i   Vectors  from  Points 
A  &  B 
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X 


Figure  A5 .   X-Z  Plane 


Figure  A4 


Y 


Plane  Formed  by  B  Vector  and  Y  Axis 

G  • 


Figure  A6 .   Plane  Formed  by  A  Vector  and  Y  Axis 


fi2 


the  inertial  cartesian  coordinate  system.   Figure  A3  is  the 
superposition  of  the  i   unit  vector  from  points  A  and  B  of 
Fig.  A2 .   The  intent  here  is  to  express  the  change  in  the 
unit  vector  i   in  terms  of  components  along  the  original 
(s,n,m)  directions  shown  at  point  A  in  Fig.  A3.   Vector  OA 

in  Fig.  A3  is  the  original  i   direction  and  vector  OB  is 

8i  ^ 

i   + ds.   Points  C  and  E  are  the  projection  of  points  B 

s    (3  s 

and  A  in  the  (x,z)  plane  respectively.   Point  D  is  the  projec- 
tion of  point  C  onto  the  line  OE . 

Figure  A4  is  the  plane  formed  by  the  OB  vector  and  the 
Y  axis.   The  length  of  OB  is  unity   by  definition.   The  length 
of  BC  is  sin (G+de) . 


sin(fi+d9)   =   sin  0  cos  dfl  +  cos  6  sin  d( 


which,  since  dO  is  a  small  angle,  is 


sin(0+de)   =   sin  fl  +  cos  6  de  (a6 


The  length  of  OC  is  cos(fl+de). 


cos(9+d9)   =   cos  9  cos  d9  -  sin  9  sin  d( 


which,  since  d9  is  a  small  angle,  is 


cos(9+d9)   =   cos  9  -  sin  9  dO  (A7 
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Figure  A5  is  the  (x,z)  plane  in  which  the  angle  (f)  is 
defined.   The  i   direction  is  shown  at  point  C.   The  length 
of  CD  is  OC  sind(|)  which,  with  Eq.  (A7)  is 

[cos  9  -  sin  9  d9]sin  dcj) 
which,  since  d(p    is  a  small  ang].e,  is 

cos  G  d(|)  -  sin  9  d9  d(p 

Dropping  the  second  order  term 

CD  =  cos  9  d(j)  (A8) 

The  length  of  OD  is  OC  cos  d(J)  which,  since  d'p    is  a  small  angle, 

is  OC. 

Figure  A6  is  the  plane  formed  by  the  OA  vector  and 

the  y  axis.   CD  is  the  projection  of  BC  and  since  both  are 

vertical  lines,  BC  and  GD  have  the  same  length,  which  has 

been  determined  in  Eq .  (A6) .   The  i   and  i   directions  are 

s       n 

shown  at  point  A.   The  lengths  to  be  determined  are  AF  and  GF, 
since  these  are  the  components  of  AB  in  the  s  and  n  directions 
For  sma] 1  angles,  OD  was  shown  to  be  equal  to  OC  and  GD  =  BC , 
so  that  the  angle  GOD  is  equal  to  angle  BOC  =  G+dB.   Angle  GOF 
is  therefore  dO,  the  length  of  OG  is  unity   and.  the  length 
of  GF  is  sin  d9,  which  for  small  angles  is  d9  ,-  and  it  is  in 
the  i   direction.   The  length  of  AF  is  OA  -OF.   OA  is  unity 
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by  definition  and  OF  is  cos  d0 ,  therefore  AF ,  the  component  of 

AB  in  the  i    direction  is  equal  to  1  -cos  d0,  which,  since 

d9  is  a  small  angle,  is  zero. 

The  component  of  AB  in  the  i   direction  is  GF  or 

n 


de  i  (A9 

n 


The  remaining  component  of  AB  in  the  i   direction  is 
CD  or 


cos  e  d(})  i  (AlO) 

^   m 


Therefore,  the  vector  AB  can  be  written  as 

O  1  /\  /\ 

ds   =   d6  1   +  cos  9  dd)  i 


3s  m  ^  "  m 


^r—  ds  1   +  COS  9  ^   ds  1 
9s      n  9s      m 


and,  dividing  by  ds,  finally 

^^s      99  :   ,  _..  .  91  ^ 


1   +  cos  9  Tr-*-  1  (All) 


s       9s   n         9s   m 


Note  that  Fig.  A3  can  be  considered  also  to  describe 

a  chanae  in  the  unit  vector  i   at  a  fixed  location,  with 

s 

respect  to  time.   In  this  case,  the  vector  AB  is  equal  to 

91 

j^   dt  and  with  Eq.  (A9)  and  Eq.  (A.IO) 
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8i 

< 


IT-   ±      +   cos  6^1 
dt   n  dt   m 


(A12) 


With  Eq.  (All) ,  Eq.  (A5)  becomes 


(V-V)V   = 


3q  :   ,   2  96  ?   ,   2 


91  i 
8s   m 


(A13) 


D.   CONSERVATION  OF  MASS 

In  the  quasi-one-dimensional  differential  stream  tube, 
shown  in  Figure  A7 ,  where  i   is  always  in  the  direction  of 
V,  p,  q,  and  the  cross  sectional  area  A  are  known  at  the  center 
of  the  element  and  p  and  q  are  assumed  constant  on  any  given 
cross  section. 


dp   ds 

P^"  8s  2 

P»^[^A_______-^ 

,  8p  ds 
P-^3s  2 

8q  ds 
"^  8s  2 

A 

A 

,  8q  ds 
a  H — —  — 
^  8s  2 

8A  ds 
^  8s  2 

1 

/ 

.  —  _ 

""  "^"J 

8Ads 
^^8s  2 

■  "ds 

Figure  A7 .   Differential  Stream  Tube  of  Variable  Cross 
Section 


The  statement  of  continuity  is 


The  change  in  the  mass 
within  the  control  volume 
with  respect  to  time 


The  net  influx  of 
mass  across  the 
control  surface 
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since  the  volume  is  constant,  the  left  side  becomes 


A  ds  ^ 


For  the  stream  tube,  mass  enters  and  leaves  only  from  the 
ends  so  that  the  right  hand  side  can  be  written  as 


r    8p  ds,  r    3q  ds,  r-   9A  ds.    r    9p  ds-,  .    9q  dSi  r^       3A  ds 


which  on  expansion  becomes 


-  p  q  3^  ds  -  qA  ^  ds  -  p  A  ^  ds  -  g|  3I  ^  ds 


On  dropping  the  higher  order  terms  and  combining  with 
the  left  hand  side, 

It  -^  P^A  9?  ^  ^  I¥  ■"  Pai  =  °  ^^^^) 

Equation  (A14)  is  the  form  of  the  continuity  equation 

which  is  required  to  model  a  flow  as  being  one-dimensional 

1  9A 
with  area  change.   In  general,  however,  —  ^^ —  must  be  expressed 

^        ^  A  9s  ^ 

in  terms  of  q,  9  and  cfi . 

With  reference  to  Fig.  A8 ,  the  change  in  the  cross- 
sectional  area  of  the  stream  tube  can  be  written  as 


7c—  ds      =       [dm   +    ^s dsj  [dn    +    ^^ ds]     -    A  (,A15 

o  S  d  S  d  S 
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3dm 


n 


dn 


Figure  A8 .   Cross  Section  of  a  Diverging  Differential 
Stream  Tube 


After  expansion,  dividing  by  A  and  dropping  the  higher  order 
terms,  En.  (A15)  becomes 


1  3A   ^   J^  9dn  ^  J^  3dm 
A  3s      dn  3s     dm  3s 


(A16 


From  Fig.  A9 , 


_1_  3dn 
dn  3s 


ds   =  -^ 


l_ 

dn 


ds 


cos (^—  dn. 
dn 


-sin  (- —  dn) 

3n 


which,  for  small  angles,  becomes 


J^  3dn   ^  dB_ 
dn  3s       3n 


(A17) 


From  Fig.  AlO, 


1  ddm  , 

:r~  ^ —  ds 
dm.  3  s 


dm 


:a18) 
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n 


Figure  A9 .   (s,n)  Plane  of  Diverging  Differential  Stream  Tube 


»  1. 


Figure  AlO.   (s,m)  Plane  of  Diverging  Differential  Stream 
Tube 


>     X 


Figure  All 


(s,m)  Plane  and  its  Projection  onto  the 
(x,z)  Plane 
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with  reference  to  Fig.  All,  AB  =  CD,  and 


\ 


CD   =   OD  sin  [|^  dm]  (A19 

drn 


OC 

OD   =   ^ (A20 

COS  r^  dm] 
"  dm 


OC   =   ds  cos  e  (A2i; 


Combining  Eqs .  (A18)  tlirough  (A21)  and  knowing  that  -r^   dm 
is  a  small  angle. 


^  |dni   .   cos  e  1^  (A22) 

dm  ds  8m  ■ 

Equation  (A16),  with  Eqs.  (A17)  and  (A22),-  becomes 


T  ^^—      =      ^T-   +    cos  e  7r±  A2  3 

A  9s      9n  9m 


With  Eq .  (A23) ,  the  general  form  of  the  continuity  equation 
in  natural  coordinates  is 


9p 
9 


^  -  <5  If  ^   p|f  ^  pqlff  -  COS  e  ||]   =   0  (A24 


E.   CONSERVATION  OF  MOMENTUM 

The  statement  of  conservation  for  an  arbitrary  control 
volume  which  is  fixed  in  the  reference  frame  can  be  written 
as 
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The  change  in  the 
momentum  of  the  fluid 
in  the  control  volume 
with  respect  to  time 


The  net  influx 
of  momentum 
across  the  con- 
trol surface 


+ 


The  net 

body 

force 
on  the  fluid 


The  net  force 
on  the  con- 
trol surface 


Assuming    the    fluid    is    inviscid,    the    only    forces    acting   on    the 
control    surface    are    pressure    forces.       In    the    absence    of    any 
body    forces     (gravity,    electromagnetic,    etc.),    the    statement 
becomes 


[^   ///    pV  dvol      =      -    //    VpV'da    -    //    P   d^ 


=   0 


(A25) 


With  Gauss'  Theorem,  Eq .  (A25)  becomes 


///  9^[pV]dvol  +  ///V- [VpV]dvol  +  ///  VP  dvol   =   0 


or 


///{^[pV]    +    V-  [VpV]    +    VPldvol      =      0 


(A26 


Since  Eq.  (A26)  is  true  for  any  volume,  the  integrand  must 
be  zero,  thus 


^[pV]  +  V-  [VpVJ  +  VP 


(A27 
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Expanding  the  first  term  of  Eq.  (A27)  and  using  the  vector 
identity 

V' [VpV]        =       V[V-pV]     +    p(V-V)V 

Eq.  (17)  becomes 


p  1^  +  V{|^  +  V-pV}  +  P(V-V)V  +  VP   =   0  (A28 


The  second  term  of  Eq.  (A28)  is  zero  by  the  conservation  of 
mass  so  that  Eq .  (A28)  reduces  to 


P  1^  +  p(V'V)V  +  VP   =   0  (A29) 


Using  Eqs .  (A2),  (A5),  (All)  and  (Al2),  collecting  terms  and 
equating  each  vector  component  to  zero,-  Eq.  (A29)  becomes 


^s[p|f  -    P^lf   -Hi       =      0  <'^30, 


^    r       99  ^     ^2   ae  ^   3p,  .  /^oi\ 

i„[pq    -at    ^    Pg       -9i    +    -9^^        =       0  (^31) 


i^[    pqcose    ^    +    pq      coseg^    +    ^]        =       0  (A32) 


Using    the    perfect   gas    relation 

a2    =    il 
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and   the    identity 


3P 

9s 


=      P 


9   enP 
9s 


Eqs .     (A30)    through    (A32)    become 


9t 


9s 


A       9  £nP 
Y     9s 


(A33) 
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96 


9t    ■"   ^    97      = 


A       9  .enP 
yq    9n 


:a34) 


9t 


A' 


9  £nP 


9s 


yq  cos  9    9m 


(A35; 


which   are    in    the   desired    form. 

F.       CONSERVATION    OF    ENERGY 

The    conservation   of    energy    for    an    arbitrary    control    volume 
is 


The  rate  of  increase  of 
energy  within  the  con- 
trol volume  with  respect 
to  time 


The  rate  at  which  energy 
is  entering  the  control 
volume  across  the 
boundary 


The  net  rate  of 
work  done  on  the 
fluid  at  the 
boundary 


:a36: 


Neglecting  gravitational  potential,  the  energy  per  unit  mass 
is 
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2 


where  e  is  the  internal  energy  per  unit  mass  and  q  is  the 
velocity  magnitude.   For  an  inviscid  fluid  in  the  absence  of 
body  forces,  the  only  work  done  on  the  fluid  is  pressure  work, 
Thus  Eq .  (A36)  is  written  mathematically  as 


2  2   _  _ 

j^   jj j[e  +^]p   dvol      =      -//[e  +^1  pV-d^  -  //  PV-da         (A37) 


Using  Gauss'  Theorem,  Eq .  (A37)  is  written  as 
2  2 


j j jj^ie  +^]p   dvol      =      -    ///v  [  (e  +^pV]dvol 

-    ///    V. [PVjdvol 


or 


2  2 

|^{  [e  +^]p}    +    V-{[e  +^]pv}      =      -^-[PV] 


Expanding , 


2  2  _  2  a^  - 

~[e  +  Y^     +    p    |^[e  +^]    +  (PV)  V-  [e  +  ^]     +     [e  +-^]  [V-  (pV)  ] 


=       -V- [PV] 


or 
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2  2  2 

[e  +^]  [|f +V-  (pV)]  +  p  |^[e  +^]  +  (pV)V-  [e  +  ^] 


=  -V- [PV]  (A38) 

The  first  term  of  Eq.  (A38)  is  zero  according  to  the  conserva- 
tion of  mass.   Expanding  the  remaining  terms  of  Eq .  (A38) , 


2  2 


or 


2 
§^[e  +2^]   =   -  i  V.  [PV]  (A39) 


Equation  (A39)  can  be  written  as 


2 


De  ^  2_[^]        1  '-  — >    P 


Dt  ■  Dt-  2^        p(V-"^)  -  p("-^)  (^^°) 


With  V  =  qi   in  the  natural  coordinate  system,.  Eq .  (A29) 
becomes 


or 


9V  ^  3V  1  ^^ 
^T^r+qTT—  =  -  —  VP 
8t     -^  9s         p 


^     =      -  i  VP  (A41 

Dt         p 
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Taking  the  dot  product  of  V  with  Eq.  (A41) , 


V§|   =   -  ^[V.VP]  (A42 


and  substituting  Eq.  (A42)  into  Eq=   (A40) 


Dt    Dt-  2  '         Dt    p^   ^^ 


But 


^-              D(qi  )      T^   -       Di 
DV   _      -^  s    =   £3.  i   +    S 

Dt        Dt        Dt   s    ^  Dt 


^   ^       9i        8i 
^q  -;    ,    r    2;   ,        St 

=   Dt  's  -^  ^f^t-  -^  ^  J^^ 


and  using  Eq  -  (A.ll)  and  Eq .  fAl2), 


Dt   =   ^  I3  +  qf^  1^  +  cos  6  3^  1^] 


+  q  [r—  1   +  cos  9^1 
^   9s   n  9s   m 


Therefore 


(A43) 


•   /Dq  ?   ^  r   99  ^   2  99.:   _,       ^  rM  ^   Mi  ■  1 
^   ^\s-^Dt  S  ^  t^  9t  ^  ^   9¥^^.  -^  ^^°^Qf9t^^  gf^^m^ 


-^  Dt       Dt   2 


Therefore  Eq .  ( A4 3 )  reduces  to 
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De        P    — 

^   =   -  ^(VV)  (A44) 


From  the  conservation  of  mass 


1^  +  V- (pV)   =   0 


or 


1^  +  V- Vp  +  p(V-V)   =   0 


which,  in  the  natural  coordinate  system  gives, 


If  .qff   =   -P(V.V, 


or 


i§^   =   (V-V)  (A45 


Substituting  Eq .  (10)  into  En.  (9)  , 


De  _  ^  Dp   ^ 

Dt     2  Dt      ^  ^^^^' 

P 


The  definition  of  the  modified  entropy  is 


dS   =   -  —  — =— 
Y   T 
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and  the  first  law  of  thermodynamics  for  an  elemental  mass 
requires  that 


dQ^   =   de  +  Pd(^) 
R  p 


so  that 


dS   =   -  ^[de  +  Pd(^)] 


:^[de  +  ^  dp] 


and  Eq.  (A46)  implies  that 


§1   =   0  (A47) 


which  states  that  modified  entropy  is  conserved  along  streamlines 

G.   MODIFIED  ENTROPY  EVALUATION 

With  the  modified  entropy  defined  as 

^Qr 

dS   = ^ 

yt 

or 


1     r^     ^^R 

S3  -S^   =   -  i   /   ^  (A48) 


From  the  first  law  of  thermodynamics 
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dQ^   =   de  +  P  dv  (A49 


where  e  is  the  internal  energy  and  v  is  the  specific  volume 
Substituting  Eq .  (A49)  into  Eq ,  (A48) , 


R  R  R 

c    c   -   1   r   de  -i-P  dv  _   1,  r   de     ,      P  dv,       ,  -. 


For  a  perfect  gas 


P  V   =   R^T       and       de   =   C   dT 

Vj  V 


for  which  Eq .  (A50)  becomes 


■   T    B  C   dT      B  R„  dv 
^A    B      Y^^^     T      ^^     V   J 


which  after  integration  yields 


T  V 

Stv   =   So  +  -[C   iln  -^  +  R  £n  — ]  (A5i: 

A       B    Y   V     T^         v^^ 


Substituting 


Rq   =   '=v''-l> 


into  Eq ,  (A51) 
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Pi  Pi  /\ 


and    since    R        =    R 


Sa      =      S3.^^^[.n!^-    y£n^] 


or 


S  S„  ,  P„  P, 

A  B    ^  1         r  „  B  o„       A    T 


or 


^B  ^\  1  ^A  ^B 

^      =      ^   +       ,       T  ,  [Jin — ^   -    in   ^-]  (A52) 

^G  ^G         y^y-'-^  O^y  PgY 


If  the  reference  condition  is  chosen  to  be 


'^   -    2        1    ,n  '* 


Rq     y-1   y(Y-i)     Pj^Y 


Eq .  {A52)  becomes 


^   =   ^ r-^^  5.n  -?-  (A53) 

R(^      Y-1    Y(Y-1)      PgY 


80 


H.   FURTHER  TRANSFORMATIONS 

Equations  (A24)  ,  (A33)  through  (A35)  and  (A47)  are  in  the 
required  form.   However,  some  of  the  variables  need  to  be 
transformed,  namely,  derivatives  of  9.n   P  need  to  be  expressed 
in  terms  of  A  and  S  and  derivatives  of  p  must  be  expressed  in 
terms  of  q  and  A. 

Using  the  first  law  of  thermodynamics  expressed  as 


dQ   =   dh  -  -  dP 

P 


Eq.  (A4  8)  may  be  written  as 


-Y  T  VS   =   Vh VP 

P 


or  for  the  i   component 


rtC       T 

9S  1         p  P        B£n    P 


8s  T       dS  pT       as 


(A5  4 


Assuming  C   to  be  constant  and  usina  the  equation  of  state 
for  a  perfect  gas  Eq.  (A54)  becomes 


3S   ^   ^  9  .  P  .  _    9  £n  P 
^  3s      T  as^R^p-"    ^   3s 


and  since  C   =  yR/(y-1 


9S_   ^   ^^G    P^G  ^r^^n  _  p   3  Zn    P 
~''  9s      (y-1  )   P   9s  R^p      G   9s 
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since  R  is  a  constant  it  may  be  moved  in  or  out  of  the 
derivatives  giving 


I     -ir-  L  T-,    J  ,,_n    -p,    .>_  L    J 


as'R  '      y-l  P  ds'p'  3s 


and  since  y  is  also  constant,  similarly 


^  9_r_S_i   =   _J P.  i_rZli    9  ^n  P 

^  9s  ^R^-"      y-l  Py  8s   p  -■  ~   8s 


.2 
Using  A   =  yP/p 


8_.^,   ^   ^^; ]^  9A^  _  9  £n  P 

"^  9s  ^R^      y-l  ^2  3s       3s 


2Ay     3A  _  3  9.n    P 

,   T  ,  . 2  3s      3s 
(y-l) A 


and,-  on  rearranging 


3  l-n    P   ^    2y  1  ^A      3   S^ 
3s        y-l  A  3s    ^  3s ^R- 


Substituting  Eq .  (A5  5)  into  Eq.  (A3 3)  yields 


and  the  derivative  of  In   P  has  been  removed 


(A55 


9q  ^    9q  ^  2A  3A  ^  ,2  3  .S.      ^  ,.^^. 

3t  -^  ^  3i  +  pi  3^  +  A   — [-]   =   0  (A56) 
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since 

p 

logarithmic  differentiation  yields 

2dA^dP_^d/   _dp 
A       P   /y       P 

For  a  pure  substance 

P   =   P  (  P  ,  S ) 
and  an  isentropic  process 

P   =   P  (  P )  3 
Differentiating 


dp  S      dp       p 


therefore 


Substituting  Eq.  (A58)  into  Eq .  {A57)  gives 


„  dA_  _     dp  _  dp   _    _,  dp_ 
A~^P     p~^"'p 


or 
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(A57) 


^   =    Y^  (A58) 

P         p 


J  2p   dA  ,-^rr  \ 


Along  a  streamline 


dp    =    It  +  g  I?  (^60) 


and  since  the  fluid  is  the  same  from  streamline  to  streamline 
Eqs .  (A59)  and  Eq .  (A60)  can  be  combined  giving 


9_P     i£  =   ^P_  dA   ^      2p  rdA  d_A. 

dt        ^    ds  Y-1  A      (Y-DA^t  ^  ^  ds- 


Substituting  Eq .  (A61)  into  Eq.  (A24)  yields 


3A  ^    9A  ^  y-1   3q      ,  Y"!  r^Q  ^     n  Mi 


and  the  derivatives  of  p  have  been  removed. 

The  governing  equations  are  now  summarized 


(A61) 


(A62) 


^-   +    q   ^-   +    -^5—  At-^   +    q  A-^5—  ^   +    cos    6    ^        =       0  (A63) 

dt  ds  2         cis  2       dn  dm 


|£  4-    q   1^  +    ^   |A   +    A^    ^[|]       =       0                  '       ■                                   (A64) 

9t  "    3s  Y~l    9s                 9s    R 

96^99  A^9PnP  ,^c^^ 

^rr  +   q  T~     =      " ^ (A65) 

9t  ^    9s  Yq         9r\ 

^"P  ^        11  A^         9    &n    P                                                                     ,-... 

^  +    q    ^      -       -    K   ^ (A66) 

9t  ^    dS  Yq  cos  9          9m 
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II  +  q  II   =   0  (A67) 


Note  that  Eqs .  (A65)  through  (A67)  are  in  the  form 


at  ^  ^  8¥  -  ^ 


If  X  =  X(s,t) ,  then 


dX   =   II  dt  +  1^  ds 

dt  dS 


and 


dx  _  dx      ds  ax 

dt   -   at  ^  dt  9¥  ^^^^^ 


where  (ds/dt)  is  the  "characteristic"  direction  in  the  (s,t) 
plane  along  which  Eq ,  (A68)  holds.   Along  this  direction ,  the 
equations  may  be  transformed  from  partial  to  ordinary  differ- 
ential equations.   V7hat  is  desirable  is  to  find  transformed 
variables  in  terms  of  which  Eqs.  (A63)  and  (A64)  will  be  in 
the  same  characteristic  form.   For  this  reason  the  modified 
Riemann  variables 


Q   =   q  +  AS/Pp  (A69 


R   =   q  -  AS/R^  (A70) 
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are  defined,  where  S/R„  is  defined  in  Eq .  (A53) .   The  modified 
Riemann  variables  can  be  Introduced  by  manipulation  of  the 
equations  as  follows.   First,  multiplying  Eq.  (A63)  by  S, 
gives 

s|4-^s|A.l^As||.,Asl^l|l.cosfi|l,  =  0 

(A71) 
and  multiplying  Eq =  (A67)  by  A  gives 


*  II  +  1  *ll  =   °  f'^"' 


On  adding  Eqs .  (A71)  and  (A72)  to  Eq,  (A64),.  and  introducing 


A  If  -  A  If   =   0  (A73) 

9s       dS 


and 


A  S  14  -  A  S  14   =   0  (A74) 

dS  dS 


after  appropriate  rearrangement 


9q    A  9S    ^  9A      da  ^  dS  ^  9A    .  3q 


9t  '  ^  3t 


+  S^^+q^+qA^-+qST—  +A^ 
c>t    -'9s    -'9s    -'9s      9  s 


^A^S-^sfA      =      -lZiAs||.As|A.A|| 


2A     9A         Y-1         ,  o  r90    -.  a    iii 

— r  ^^  ~   ^r-  q  A  S  hr-   +    cos    6    -;r±        =      0 
■y-l    9  s  2       ^  9n  9m 
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Rearranging  and  introducing  Eq.  (A69)  gives 


9Q  ^  ,    ^A^     3Q        7-1  Arc    2  T  r 9  /     2  ,. , 


Yzlg,S[f  .cos    e|l]  (A75) 


By  subtracting  Eqs .  (A71)  and  (A72)  from  Eq.  (A64),  again  using 
Eq.  (A73)  and  Eq .  {A74),  and  introducing  Eq .  (A70),  one  also 
obtains 


_+     (q-A)     —       =       Y-A[S    -^[^(q+^33- 


+   ^  AS[|^  +    cos    e    1^]  (A76) 

z  dn  dm 


Together  with  Eqs.  (A65)  through  (A67) ,  Eqs.  (A75)  and  (A76) 
form  the  system  of  equations  which  describe  the  isentropic  flow 
of  an  inviscid  perfect  gas  under  unsteady,  compressible  condi- 
tions, and  which  may  be  solved  as  a  system  of  ordinary  differ- 
ential equations  along  characteristic  trajectories  in  the 
space-time  domai.n. 

I.   SHOCK  JUMP  EQUATION 

The  analytical  expression  for  the  change  in  the  extended 

Riemann  variable,  Q,  across  a  stationary  shock  is  derived 

below. 

Consider  a  shock  movina  with  a  velocity  V   into  a  fluid 

^       s 

with  velocity  qj,  as  dep.i  cted  in  Fig.  A12,  with  the  conditions 


A  O- 


V 


-O  B 


u,  = 


q  -V 
^a  s 

A  O- 


-o   B 


%% 


Figure  A12 


Figure  A13 


to  the  left  of  the  shock  being  denoted  by  A  subscripts  and 
conditions  to  the  right  by  B  subscripts.   If  a  velocity  of  -V 
is  imposed  on  the  entire  system,  the  situation  becomes  one  of 
a  stationary  shock,  as  depicted  in  Fig.  A13.   In  both  cases, 
the  high  pressure  side  is  to  the  left  (A  side) .   Since  all 
velocities  are  defined  positive  to  the  right,  a  positive  value 
for  the  relative  mach  number  is  defined  as 


u. 


W   =   - 


B 


(A77) 


The  extended  Riemann  variable,  Q,  is  defined  as 


q  +  As 


:a78) 


where 


-1 


Y  (Y 


hr  'nt^, 


(A79) 


If  all  velocities  are  non-dimensionalized  by  A^^  and  pressures 
and  densities  are  non-dimensionalized  by  their  values  at  B, 


^A  -  Qb  ^A-^B    ^    ^    ,  , 


^A-^B  ^  ^A.   2  1       „    Pa  Pa/ 


2  1  .     Pb.^B.Y 


5-n-^(-^)     ] 


Y-1      y(y-1)     '     Pg    Pg 


"a      %    ,    /a       ,.     2  /a,         1         ro       ^A /b.y, 

=         ^ +    ll 1]  T     "       (-1—)  } ;-^L  ^n     (— r-)        ] 

Ag^        Ag         'y-1  Ag    Y(Y-l)  Pg    pA' 

(ABO) 

The  ratios  of  pressure,  density  and  sound  speed  across  a 
normal  shock  are  well  known  functions  of  Y  and  the  Mach  number 
which  in  this  case  is  W.   They  are 


A       2y   2    Y-1 

— tr   Ta; 


Pg      Y  +  1  "^  Y  +  1 


^A   _     (y  +  1)  W^ 


Pb      (y-I)W^  +2 


;a82 


b 


An  expression  for  (q^  -q^)/A   is  obtained  using  mass 

ABB 

conservation . 
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A  A      ^B  B 


"a   _   ^   _   (y-1.)v7^  +2 
^B      ^A  (Y  +  1)w 


Subtracting  one  from  both  sides. 


"^A   ^B   ^   (y-I)w^  +2  -  (Y  +  l)w^ 
""b  (Y+1)W^ 


A    B   _   2  (1  -  w  ) 


,   .  ^  2 

B         (Y+1) W 


Substituting  Eq.  {All)     for  u   in  Eq .  (A84)  gives 


^A   '^B   ^   2  ( 1  -  w^  ) 
B        (Y+1)  W 


''a  -  ''b      _       2(w^--1) 
Ag        (y+1)  W 


Since  u,  =  q,  -V   and  u_  =  Q^^  -  V 
A    ^A    s       B    ^B    s 


^A-^B   =   ^A-^s  -  ^B  ^  ^s   =   ^A-^B 


Thus  Eq .  (AB5)  can  be  written  as 


(A84) 


(A85) 
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'^A-% 
~^^ 


^      2(w  -1) 
(Y+1)W 


(A86) 


Substituting  Eqs ,  (A86)  and  (A81)  through  (A83)  into  Eq .  (A80) 
gives 


^a-Qb 


2 

W-1)_L  2     r        1    n,   Tx/i,Y~l-2..2y2 


2_[W   -1) 
(Y 


1/2 


-  1}  -  {t-^[2(y-i)(i+^w')(-|^w'-1)]'/^> 


2  ^  Y 


X  (-l^Un(^w'-^)(^^-^^^  ;^)'} 


Y  (y-1)     Y  +  1     Y  +  1 


;a87) 


(y+i)w' 
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APPENDIX  B 
USING  EULER-1  ON  THE  NPS  VM/CMS  SYSTEM 

A.  MEMORY  REQUIREMENTS 

Storage  should  be  defined  as  one  mega-byte  to  ensure  ade- 
quate memory  for  DISSPLA  requirements.   Also  ensure  that  room 
exists  on  the  user's  disk  for  the  listing  file  if  that 
output  option  is  selected. 

B.  TERMINAL  REQUIREMENTS 

The  type  of  terminal  required  to  run  EULER-1  depends  on 
the  output  option  selected.   When  using  the  tabular  listing 
output,  any  terminal  tied  into  the  VM/CMS  system  may  be 
used.   Graphical  output  may  also  be  created  using  any  terminal 
and  the  graph  stored  on  the  user's  disk  in  a  metafile  for  later 
viewing  at  a  graphics  terminal  or  for  printing.   To  have  the 
graph  stored  on  the  disk,,  use  the  "COMPRS"  command  on  line 
223  of  the  program,  and  comment  out  the  "TEK618"  command  on 
line  224.   When  graphical  output  is  selected  and  it  is  desired 
to  view  the  graph  at  the  terminal  at  execution  time,  an  IBM 
3277-Tek618  dual  screen  terminal  must  be  used.   To  use  this 
option,  use  the  "TEK618"  command  on  line  224  and  comment  out 
the  "COMPRS"  command  on  line  223.   VJhen  using  this  option,  a 
low  quality  hard  copy  of  the  TEK618  display  may  be  obtained 
using  the  TEK4631  hard  copy  unit  attached. 
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C.   PROCEDURE  AND  COMMANDS 

After  loging  on  to  VM/CMS ,  use  the  command  "DEFINE 
STORAGE  IM"  to  increase  the  virtual  memory  as  recommended 
above.   This  will  remove  you  from  the  CMS  environment.   Re- 
enter CMS  with  the  command  "I  CMS." 

Edit  the  user  input  area  of  EULER-1  FORTRAN  as  desired 
for  the  particular  problem  being  run.   When  graphical  output 
is  selected,  comment  out  either  the  "COMPRS"  or  the  "TEK618" 
comjn:\and  as  desired. 

Compute  the  program  using  the  command  "FORTVS  EULER-1." 
After  compilation,  an  EULER-1  TEXT  and  an  EULER-1  LISTING 
file  will  reside  on  the  user's  disk.   At  execution,  if  a 
tabular  output  has  been  selected,  it  will  be  sent  to  the 
user's  disk  as  "FILE  FT09F001,"   To  give  this  file  a  particu- 
lar name,  say  "EULER-1  LISTING,"  use  the  file  definition 
command 

"F  I  9  DISK  EULER-1  LISTING  A  (  PERM" 

at  compile  time.   An  alternate  and  convenient  means  of  com- 
piling the  program  is  to  set  up  an  exec  file  on  the  user's 
disk  by  the  name  "EULER-1  EXEC"  which  contains  the  commands 

F  I  9  DISK  EULER-1  LISTING  A  (  PERM 
FORTVS  EULER-1 

Then  to  compj.le  the  program  and  define  the  output  file,  use 
the  command  "EULER-1." 

After  compilation,  to  execute  the  program,  use  the 
command  "DTSSPLA  EULFR-1."   The  message 
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...YOUR  FORTRAN  PROGRAM  IS  NOW  BEING  LOADED... 
...EXECUTION  WILL  SOON  FOLLOW... 

should  appear,  followed  shortly  by 

. . .EXECUTION  BEGINS. . . 

If  tabular  output  was  selected,  proper  termination  will 
result  in  a  ready  message.   If  graphical  output  was  selected 
using  the  TEK618  option,  the  plot  will  appear  on  the  TEK618 
terminal,  the  program  will  pause,  and  a 

"...press  ENTER  to  continue" 

message  will  be  displayed  on  the  3277  terminal.   At  this 
point  the  hard  copy  must  be  made,  if  desired,  using  the 
TEK4631  hard  copy  unit.   After  pressing  the  ENTER  key  on  the 
3277  terminal,  the  plot  will  be  erased  and  the  program  will 
terminate.   Proper  termination  will  be  indicated  by  the  m.essage 

"END  OF  DISSPLAY  9.2    ####  VECTORS  IN  1  PLOT,.." 

and  a  ready  message. 
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APPENDIX  C 
EULER-1  FORTRAN  CODE 

OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO 

0  0  000  00  00•-^-^'-H--4-^•-^--lr-^-^,^(^Jruf^l^^J^^JfN^^^vJC^Jf^J^O^O^o^O'^^Oc^rn^O'^vt  Nj-<f"  <f  <i- -4-  ^t  ~t  vj- 

oooooooooooooooooooooooooooooooooooooooooooooooo 
oooooooooooooooooooooooooooooooooooooooooooooooo 

^_J_I_J_J_J_J_J_J_J_J_U_I— I^_J_I_I_J  J_U^-J_J  J-J— l_l-J_i— l_)U_J-J  J— l-J_J_l_i-J_J_J_l_l 

r^ -^ -5  ^5 -5  ")  r) -:3  T3  rD -^  "3  r) -5  "5 -) -5  r> -5  "5 -5 -)  rD -n -5 -5 :::) --^  "^ -5  "r^ -i  ^3 -3  rj -3  ri -i  z) -) -)  ::5 -D -5 -5  r5 -)  r) 
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1 
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o> 

* 

OrO 
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1 
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* 
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